Magnetic Phases of Two-Component Ultracold Bosons in an Optical Lattice 
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We investigate spin-order of ultracold bosons in an optical lattice by means of Dynamical Mean- 
Field Theory. A rich pheise diagram with anisotropic magnetic order is found, both for the ground 
state and at finite temperatures. Within the Mott insulator, a ferromagnetic to antiferromagnetic 
transition can be tuned using a spin-dependent optical lattice. In addition we find a supersolid 
phase, in which superfiuidity coexists with antiferromagnetic spin order. We present detailed phase 
diagrams at finite temperature for the experimentally realized heteronuclear ^'^Rb - mixture in 
a three-dimensional optical lattice. 
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I. INTRODUCTION 



Ultracold atoms in optical lattices give access to stud- 
ies of quantum magnetism with unprecedented control 
and precision. Whereas efficient cooling of fcrmionic 
atoms in optical lattices remains an experimental chal- 
lenge, mixtures of bosonic atoms can more easily be 
cooled to the relevant temperature scales, thus realizing 
the Bose-Hubbard ModeUii^ with multiple species. In- 
deed, recent experiments have already succeeded in load- 
ing a heteronuclear mixture into an optical lattice^. The 
interaction between the two species has been adresscd 
by means of a Feshbach resonance^, which offers a di- 
rect way of mapping out the phase diagram as a function 
of the interspecies interaction strength. Moreover, su- 
perexchange processes have been directly observed in a 
mixture of two-component bosons^. 

Whereas the phase diagram of spinless bosons is 
qualitatively well captured by Gutzwiller Mean-Field 
theory^, for a multispecies system dynamical correla- 
tions are important, because they give rise to spin or- 
der. Previous theoretical studies have either discussed 
the weak tunneling limitii^ or performed expansions 
around mean-field^ or the strong-coupling limi t^^'^^ . 
Numerical non-perturbative methods are available in 
one spatial dimensio n^^d^ , and Quantum Monte-Carlo 
simulations were very recently performed for two spa- 
tial dimensionsi^. In this article we focus on the 
generic three-dimensional situation where Dynamical 
Mean-Field Theory (DMFT) has been estabHshedi^ as 
a highly reliable, nonperturbative approach to strongly 
correlated fermionic quantum systems. Here we apply 
the bosonic version of DMFT (BDMFT), as recently in- 
troduced by Byczuk and Vollhardti^. BDMFT treats 
condensed and normal bosons on equal footing. It is 
non-perturbative and hence can be applied within the 
full range from small to large couplings. The control pa- 
rameter is the lattice coordination number z, i.e. the 
theory becomes exact in infinite dimensions. We present 
the BDMFT equations as a controlled 1/ z expansion, up 
to subleading order. Our derivation is thus different from 
the original proposali^, in which BMDFT is constructed 



as a well-defined theory in strictly infinite dimensions, 
which requires a different scaling of superfluid and nor- 
mal parts of the action. In contrast, our derivation is 
based on a uniform scaling ^ 1/z of the bosonic hopping 
amplitude. To leading order this yields Gutzwiller Mean- 
Field theory^, while from the subleading terms of order 
0{\/z) we obtain the BDMFT equations. We thus re- 
gard BDMFT as an expansion in 1/z around Gutzwiller, 
which in our opinion is the most natural viewpoint. 

In practice, BDMFT turns out to be an efficient and 
fast scheme, which allows to map out phase diagrams 
with high resolution. Moreover, it not only allows for 
non-perturbative calculation of local obscrvables, but 
also the spectral function, relevant for RF-spectroscopy, 
is directly accessible. 

While in Ref. [3 calculations were only performed for a 
simplified lattice model with partially immobile bosons, 
here we apply bosonic DMFT to the full two-component 
Bose-Hubbard model in finite spatial dimensions. Be- 
sides the superfluid, we identify Ay-ferromagnetic and 
Z-antiferromagnctic phases, in which translational sym- 
metry is spontaneously broken and anisotropic magnetic 
order arises. Moreover, we find a supersolid phase where 
superfluidity coexists with antiferromagnetic spin order. 
We investigate the stability of these phases against ther- 
mal fluctuations, paying special attention to the experi- 
mentally relevant case of a heteronuclear ®^Rb - ^^K mix- 
ture in a three-dimensional optical lattice^. 



II. METHOD 

We first describe how BDMFT can be implemented for 
the Bose-Hubbard Hamiltonian. In particular we iden- 
tify an Anderson Hamiltonian that reproduces the local 
effective action and allows us to solve the BDMFT self- 
consistency problem. 

The starting point for our investigation is the multi- 
species single-band Hubbard model within tight-binding 
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approximation 

H ^-Y^ (t^bl^bj^ + h.c.) +\Y1 f^^M^ii. ('V " ^''^^ ' 

. (1) 

which provides an accurate description of bosonic atoms 
in a sufRciently strong optical lattice. Here t^, are 
(species-dependent) hopping amplitudes, U^^ contains 
the inter- and intraspecies interactions and (ij) indicates 
a summation over nearest neighbor sites i and j. In the 
spirit of the cavity derivation of the fermionic DMFT 
equationsi^ we consider a single lattice site (called the 
"impurity site" ) and formally integrate out all the other 
degrees of freedom. This defines the effective action of 
the impurity site as 

^imp = 4fi)= j n^'^S.^^'^o^.e-^--, (2) 

where Z is the full partition function and Z^^^ is the par- 
tition function of the cavity system without the impurity. 
For reasons of brevity we derive the effective impurity ac- 
tion and perform the numerical calculations in this article 



for the case of a Bethe lattice (Cayley tree), which in in- 
finite dimensions (z — > oo) has a semicircular density of 
statesi^: /Oo(e) = \/4zt^ — e^/27rzi^. The use of the semi- 
circular DOS in our calculations has merely technical 
reasons because this choice simplifies the DMFT equa- 
tions. Our obtained results remain qualitatively similar 
for any symmetric DOS representing a bipartite lattice. 
This is in particular true for the three-dimensional cubic 
lattice, which has only mild Van Hove singularities. For 
fermionic DMFT it has been established that the agree- 
ment between results on the Bethe lattice and the cubic 
lattice is not only qualitative, but also quantitative, with 
a typical accuracy of around ten percent. We find that 
the same is true for BDMFT, as we show below for the 
case of single component bosons, where we compare the 
BDMFT results with numerically exact Quantum Monte 
Carlo results. 

In deriving the effective impurity action, we first for- 
mally rescale all hopping parameters as = t^jjz^ such 
that 1/z appears as the small parameter in the theory. 
Based on the linked cluster theorem, the action of the im- 
purity site up to subleading order in 1/z is then obtained 
in the standard wayi^ as: 
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Here we have defined 

K.{r) = {brAr))^ (4) 
as the superfluid order parameters, and 

GUATy)^-{hAr)blAr%^ + MrW,A^'), (5) 
G^.,,,(T,r')=-(6,,^(r)6,-,(r'))o + </'.,.(r)</>,,p(T') (6) 

as the diagonal and off-diagonal parts of the connected 
Green's functions, respectively. The notation (. . .)o 
means that the expectation value is taken in the cav- 
ity system excluding the impurity site. For finite z the 
action Q coincides with the one previously derived 
Note however, that our derivation is different. In the 
original proposali^, bosonic Dynamical Mean-Field The- 
ory (BMDFT) is constructed as a well-defined theory in 
strictly infinite dimensions, which requires different scal- 
ing of superfluid and normal parts of the action. In con- 
trast, our derivation is based on a uniform scaling ~ 1/z 
of the bosonic hopping amplitude, since we focus on fi- 
nite dimensions and and our goal is to make direct con- 
tact with the three-dimensional experimental situation. 



The terms involving Green's functions in the action ([3]) 
are of order 0(1/ z), since they come with two factors 
oi ty ~ 1/z and one summation over neighboring sites, 
which gives a factor z. All the other terms are of order 
0(1); for the last term in the action this follows from 
the fact that it involves one factor oi ~ 1/z which 
cancels against the factor z arising from the summation 
over neighboring sites. Therefore to leading order the ac- 
tion ^ yields Gutzwiller Mean-Field theory^, while by 
including the subleading terms of order 0(l/z) we ob- 
tain the BDMFT equations. Hence we regard BDMFT 
as an expansion in 1/z around Gutzwiller; this is in our 
opinion the most natural viewpoint. 

To proceed, expectation values in the cavity system 
need to be identified with those on the impurity site, in 
order to obtain a closed self-consistency loop. Since sites 
at the edge of the cavity have one neighbor less compared 
to the impurity site (see Fig. [1]) , simply identifying the 
expectation values yields an error of order 1/z. For the 
Green's functions this poses no problem, because they 
already appear at subleading order in the action, but it 
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FIG. 1: Illustration of the cavity method. Sites which are con- 
nected to the impurity (coloured greenish) have one neighbor 
less once the impurity is removed. 



yields a relevant correction to the superfluid order pa- 
rameter and turns out to be essential for quantitatively 
accurate predictions of the phase diagram. Details re- 
garding the implementation will be given below. 

We now turn to the solution of the effective action. 
This we do in the spirit of the exact diagonalization (ED) 
solution of fermionic DMFT— . We represent the effective 
action ([3]) by an Anderson Impurity Hamiltonian H^: 



Ha = zty{(l)lb^ + h.c.) + ^^Ui,f,n^{n^- 6^^,) 



E 



(7) 



The chemical potential and interaction term are di- 
rectly inherited from the Hubbard Hamiltonian. The 
Gutzwiller term represents the bath of condensed bosons 
with superfluid order parameters 4>i, for every component. 
The bath of normal bosons is modeled by a finite num- 
ber of orbitals with creation operators aj and energies ei . 
These orbitals are coupled to the impurity via normal- 
hopping amplitudes V^^i and anomalous-hopping ampli- 
tudes W^^i . The anomalous hopping terms are needed to 
generate the off-diagonal elements of the hybridization 
function. We define the following hybridization func- 
tions: 



I 



(8) 
(9) 



Integrating out the orbitals leads to the same effective 
action (3), if the following identification is made: 



(10) 



These self-consistency conditions are completed by the 
condition for the superfluid order parameter 



(11) 



The notation 



\z-l 



means that the expectation value 



is corrected for the missing neighbor on the sites adjacent 



to the impurity. Since this is a correction of order 0(1 /z) 
and 1/z is small, this correction is implemented by means 
of perturbation theory in 1/z. Equations (jlOp and (fTT|) 
thus consitute the set of BDMFT self-consistency condi- 
tions. 

The self-consistency loop is solved as follows: starting 
from an initial choice for the superfluid order parameter 
and the Anderson parameters, the Anderson Hamilto- 
nian is constructed in the Fock basis and diagonalized to 
obtain the eigenstates and eigenenergies. New superfluid 
order parameters are then obtained from (j)^, = (&i/)q~^. 
The eigenstates and energies also allow us to calculate the 
Green's functions. Subsequently, new Anderson parame- 
ters are obtained by fitting the hybridization functions 
to their corresponding Green's functions according to 
Eq. (jlOp. which is done by a conjugate gradient method. 
With this new Anderson Hamiltonian the procedure is 
iterated until convergence is reached. 

We note here that this derivation is independent of 
temperature. This implies that we cannot only deter- 
mine ground state properties, but also obtain information 
about the thermodynamics of lattice bosons, as we will 
show in the following results. Similar to fermionic DMFT 
this raises the question how BDMFT deals with situa- 
tions with broken symmetries, in which case Goldstone 
modes are present in the spectrum. Indeed, the gapless 
long wavelength excitations are absent from the DMFT 
spectrunii^. However, since in three dimensions the spec- 
tral weight of the Goldstone mode is finite and gener- 
ally small, this approximation can be justified and does 
not prohibit qualitative agreement between (B)DMFT 
results and more exact methods (if available), even in 
symmetry-broken states. 



III. RESULTS 

A. Single component bosons 

We now first apply BDMFT for the case of single 
component bosons, in which case we can compare the 
results with numerically exact Quantum Monte Carlo 
data— and strong coupling expansionsi^ii^ on the cu- 
bic lattice, and with the exact solution on the Bethe 
lattice^a. Solving the BDMFT equations for the single- 
component Bose-Hubbard model leads to an extension 
of the Mott-insulating lobes compared to the mean-field 
results (see Fig. The agreement with the exact re- 
sults on the Bethe lattice^ is very good: for the lowest 
Mott lobe the phase boundaries agree within a few per- 
cent, whereas for the higher Mott lobes the agreement is 
even better. The predicted shift on the cubic lattice is 
slightly larger-i2ii2iiiSj which is due to the different lattice 
structure. This quantitative agreement with the exact 
solution clearly shows that the applied 1/z expansion is 
a very good approximation for a three-dimensional sys- 
tem with z = 6. It is important to note that to obtain 
this result the 1/z correction of the superfluid order pa- 
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FIG. 2: Single component phase diagram for various temper- 
atures as obtained by BDMFT. In the inset: number fluctu- 
ations for density n = 1 and T = 0. The lattice coordination 
number is chosen as z = 6. 

rameter discussed in the previous section is crucial. 

Besides this quantitative agreement regarding the 
boundary of the Mott insulating lobes, it is also impor- 
tant to note that BDMFT predicts nonzero number fluc- 
tuations in the Mott insulator, as shown in the inset of 
Fig. [21 The number fluctations are clearly finite in the 
Mott state and decay as 1/U for large U . This is in con- 
trast to the Gutzwiller mean- field result, where number 
fluctuations strictly vanish within the Mott state. Since 
these finite fluctuations in the Mott state are essential for 
resolving spin order in the case of multicomponent mix- 
tures, spin order cannot be resolved withing Gutzwiller. 
In contrast, the BDMFT contains number fluctuations 
in the Mott insulating state and is able to describe spin 
order as wc describe in the following subsection. 



B. Two-component bosons 

The two-component Bosc-Hubbard model has a very 
rich phase diagram, because additional spin order exists 
in the Mott phase. Here we focus on the situation that 
the total particle density is fixed at one boson per site. 
In the strong coupling limit, i.e for t^, <^ U^^, the system 
can be mapped to a spin model, which predicts the ex- 
istence of a Z-antiferromagnet and a Xy-ferromagnetic 
phase^i^. In terms of the particle creation/annihilation 
operators 6i , &2 , all the insulating phases have the prop- 
erty that {hi) = (62) = 0. The Z-antiferromagnetic phase 
breaks the translational symmetry and is defined by the 
antiferromagnctic order parameter A^j = Itt-q,^^ — riQ,^| 
being nonzero, where fi denotes the component and a 
(a ~ —a) the sublattice. The correlator (6^62) vanishes 
in the Z-antiferromagnetic phase. The Xy-ferromagnet 
is defined by the local correlator (&IS2) being nonzero, 
and is also termed a counter-flow supcrfluid. This state 
does not break the translational symmetry and hence 
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FIG. 3: (Color online) Phase diagram of a two-component 
bosonic mixture (in an optical lattice with z = 4) at total 
density ni -f 712 = 1 for different temperatures and ratios of 
the interspecies to the intraspecies interaction. 



At weaker coupling the spin model breaks down, be- 
cause quantum fluctuations become important. A fluctu- 
ation calculation^ and strong coupling expansioniS have 
extended the spin model results to weaker coupling, 
where a transition to a superfluid phase (61), (62) 7^ 
takes place. The superfluid does not break the trans- 
lational symmetry (A^^j = 0), but shows XY-ordering: 

\{b\h)\>\{b\){k)\. 

We now investigate this system by means of BDMFT 
which allows us to study the full range from weak to 
strong coupling and effects of flnite temperature. We 
flrst study the case that the intraspecies interactions 
U = Ui = U2 are equal and much larger than the in- 
terspecies interaction U12 ^ U. We moreover vary the 
hopping amplitudes ti and t2. The condition nij + rid ~ 1 
is enforced by the choice of the chemical potential fii = 
II2 = U12/2. This has the consequence that the rela- 
tive density of the two components is changing through- 
out the phase diagram: the species with a higher hop- 
ping constant (the light species) has a slightly higher 
density in the superfluid phases of the phase diagram. 
In the insulating phases, on the other hand, the Mott 
gap protects the particle number and the relative den- 
sities are to a good approximation equal. Results are 
presented in Fig. [3] for various ratios U/Ui2- For easy 
comparison with Ref. [ij here z = 4 is chosen. In agree- 
ment with the fluctuation calculation^ we obtain a XY- 
ferromagnetic state when the hopping amplitudes are 
comparable. This XF-ferromagnetic domain shrinks if 
U /U12 becomes larger. For a larger difference between 
the hopping amplitudes there is a flrst order phase tran- 
sition towards a Z-antifcrromagnctic state. The XY- 
ferromagnetic to supcrfluid transition is of second order. 
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FIG. 4; (Color online) Phase diagram of a *^Rb-''^K-mixture (for z = 6) at fixed total density + uk = 1 as a function of 
lattice depth s and Rb-K-scattering length in units of the Bohr radius. 



while on the other hand the Z-antiferromagnet to super- 
fluid transition is first order, which is the reason for the 
large coexistence regions where both Z-antiferromagnetic 
order and superfluidity are stable/metastable, depend- 
ing on the initial conditions. Within these metastable 
regions, patches with ground state spin order will nu- 
cleate into the metastable background. Therefore this 
phase will be susceptible to the appearance of quantum 
emulsions, as predicted for a one-dimensional Bose-Bose 
mixturo^i. However, this phenomenon is expected to be 
much more prominent in one dimension than in the three- 
dimensional situation we consider. 

For very anisotropic hopping amplitudes, we predict 
the appearance of a supersolid phase. In this supersolid 
the species with the smaller hopping amplitude (for con- 
venience we take this as species 2, i.e. ti ^ is insulat- 
ing: (62) = 0, whereas the other component (i.e. species 
1) is superfluid: (61) ^ 0. Moreover the complete sys- 
tem breaks translational symmetry by developing a spin 
density wave, in which both components have alternating 
high and low on-site densities: A^^j > , for fj, = 1,2. 

Since the light species breaks the U{1) and transla- 
tional symmetries simultaneously, this phase is a true 
supersolid. The breaking of the translation symmetry is 
mediated by the presence of the heavy species, akin to 
the situation in a Bose-Fermi mixture, where supersolid 
order has been predicted as wc\&^. The supersolid has 
a flrst order transition to the superfluid and a second or- 
der phase transition to the Z-antiferromagnct. The latter 
one can be understood as the localization transition of 
the light component. The supersolid phase has not been 
predicted in earlier studies on Bose-Bose mixtures^i^^. 
Only very recently it was observed in a QMC-analysis 
for the two-dimensional casein. 

For nonzero temperatures additional quantum phases 
appear (see Fig. [3]). At low temperatures the XY- 
ferromagnet and Z-antiferromagnet in low-hopping- 
regions develop into an unordered Mott-state with (6^) = 
{b\b2) = A|^f = 0. The coexistence regions and insulator- 
to-superfluid transitions remain unaffected. For higher 
temperatures the XY-insulating phase is reduced to a 
small strip between the growing unordered phase and 
the receding superfluid. The Z-antifcrromagnctic and 
the AF-SF-cocxistcncc region diminish considerably. For 



certain parameters the counter-intuitive phenomenon of 
reentrant superfluidity takes place: the low-temperature 
antiferromagnctic phases become superfluid when the 
temperature is increased. This is the case because the su- 
perfluid is more stable against temperature fluctuations 
than the insulating Z-antiferromagnet. If the tempera- 
ture is increased in the region of stability of the super- 
solid, first the translational symmetry is restored: we 
obtain a phase in which only the light component is su- 
perfiuid ((61) ^ 0, (62) = 0), but which has no broken 
translational symmetry (A^^j = 0) in contrast to the su- 
persolid. We call this phase monofluid. Upon further in- 
creasing the temperature, also the remaining superfiuid 
order is lost. 



C. Rubidium- Potassium mixture 

Up to now, theoretical calculations were mainly per- 
formed for the symmetric parameter choice Ut = Ud- 
However, the experimentally at present most relevant 
Bose-Bose mixture consisting of ^^Rb and '*^K generally 
does not have this property^il. Here we choose the wave- 
length of the optical lattice equal to A = TSTnm, which 
yields equal dimensionlcss lattice depths s — Vq/Er for 
the two species. Er is the recoil energy and Vb is the 
strength of the optical potential, which is proportional 
to the product of laser intensity and atomic polarizabil- 
ity. The ratio of the intraspecies interaction parameters 
is then fixed according to J/Rb/C^K = WKORb/'TiRbak ~ 
0.72. The ratio of the hopping coefficients is also fixed: 
^Kh/tK ~ 'TiK/mRb ~ 0.47. This choice of the wave- 
length turns out to be sufficiently anisotropic to show 
both XF-order and antiferromagnctic order, which is not 
possible for mixtures of different hyperfine states of the 
same atom. Choosing the wavelength far red-detuned 
like in Ref. |3 on the other hand, excludes the XY-pha.se 
from the phase diagram, but makes it possible to study 
the antiferromagnet and supersolid. 

Experimentally, the ratio of intra-species interaction 
to inter-species interaction can be tuned via Feshbach- 
resonances^. Furthermore the optical lattice depth s 
can be changed to tune the ratio J/Rb/^Rb- We inves- 
tigate the resulting s-aRbK phase diagram at fixed total 
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density riRb + n-K = 1 by means of BDMFT, still us- 
ing the semicircular DOS, but taking now lattice coor- 
dination number z 6 as corresponding to the three- 
dimensional situation. Results are shown in Fig. [31 This 
mixture displays superfluid, XY-ferromagnetic and Z- 
antiferromagnetic phases and we also observe the hys- 
teresis (metastable) region between superfluid and anti- 
ferromagnet. For nonzero temperatures and high lattice 
depths the unordered Mott-state appears. At the tem- 
perature T = 2.5 X 10~^ ER^nh, {En.nb being the recoil 
energy of Rubidium), the ordered insulating states are 
reduced to a small part of parameter space, which is di- 
minished even further for higher temperatures. In order 
to compare this temperature to the temperatures reached 
in recent experiments, we consider a simple model of free 
bosons that undergo adiabatic time evolution while the 
optical lattice is ramped up. We obtain an estimate for 
the temperature at the relevant lattice depth in the Flo- 
rence group that is one order of magnetitude larger than 
the highest temperature investigate here^. Recent direct 
measurements of the temperature of a spinful bosonic 
mixture in an optical lattice at MIT have yielded temper- 
atures which correspond to only twice the highest tem- 
perature we consideri^l. 

The phase diagrams in Fig. [4] are valid for the case of 
a shallow harmonic trap, where in the trap center the 
potential is very flat. Moreover, the two species need to 
be equally distributed in the trap center, which means 
that the gravitational sag has to be compensated. 



IV. CONCLUSIONS 

We have derived bosonic DMFT within a 1/z- 
expansion and extended the formalism to the full multi- 



component Bosc-Hubbard model in finite dimensions. 
We first validated the method by applying it to spin- 
less bosons. Qualitative and quantitative agreement with 
other methods was established. Subsequently we investi- 
gated a two- component mixture, we applied the method 
to a two-component mixture. A rich phase diagram in- 
cluding spin-ordered and supersolid phases was found. 
We furthermore calculated the experimentally relevant 
phase diagrams for a *^Rb- '*^K in an optical lattice at 
zero and finite temperature. 

Note added. Recently, the BDMFT equations were also 
solved for the single-component Bose-Hubbard model in 
infinite dimensions, using a similar Anderson Hamilto- 
nian and the same Exact Diagonalization approach as 
originally proposed by us^^. 
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